#=========================================================
# Sampling Bias Analysis
# Spirostreptida millipedes in Atlantic Forest
#
#
# Purpose:
# Evaluate spatial sampling bias of millipedes of the order Spirostreptida using the package "sampbias"
#=========================================================

library(raster)
library(terra)

lat_limite_inf <- -33.8
lat_limite_sup <- 6.0
lon_limite_inf <- -56.0
lon_limite_sup <- -32.0

setwd("C:/Users/Motion/Desktop/Vieses AF dados")
resolucoes <- c(0.25) 

for (res in resolucoes) {
    r <- raster(
        xmn = lon_limite_inf, 
        xmx = lon_limite_sup, 
        ymn = lat_limite_inf, 
        ymx = lat_limite_sup, 
        resolution = res, 
        crs = "+proj=longlat +datum=WGS84"
    )
    r[] <- 1  
    nome_arquivo <- paste0("raster_", gsub("\\.", "_", res), ".tif")
    writeRaster(r, filename = nome_arquivo, format = "GTiff", overwrite = TRUE)
}

setwd("C:/Users/Motion/Desktop/Shapes/Fonte httpmapas.mma.gov.bri3geodatadownload.htm/Mata Atlântica Biomas 5000")
shape_AF <- vect("Mata_Atlantica_do_Biomas_5000.shp")


pasta <- "C:/Users/Motion/Desktop/Vieses AF dados"
rasteres <- list.files(pasta, pattern = "\\.tif$", full.names = TRUE)

getwd()
for (raster_path in rasteres) {
    r <- rast(raster_path)  # Carregar raster
    r_crop <- crop(r, shape_AF)  # Cortar raster
    r_mask <- mask(r_crop, shape_AF)  # Aplicar máscara
    nome_saida <- file.path(pasta, paste0("cortado_", basename(raster_path)))
    writeRaster(r_mask, nome_saida, overwrite = TRUE)
}

setwd("C:/Users/Motion/Desktop/Vieses AF dados")
rasteres_cortados <- list.files(pasta, pattern = "^cortado", full.names = TRUE)


for (raster_path in rasteres_cortados) {
    r <- rast(raster_path)  # Carregar raster
    r_poly <- as.polygons(r, dissolve = TRUE)  
    nome_saida <- file.path(pasta, paste0("shape_", tools::file_path_sans_ext(basename(raster_path)), ".gpkg"))
    writeVector(r_poly, nome_saida, overwrite = TRUE)
}


library(sampbias)
library(terra)
library(tidyverse)
library(data.table)
library(sf)
library(viridis)
library(purrr)
library(hablar)
library(rnaturalearthdata)
library(rnaturalearth)

setwd("C:/Users/Motion/Desktop/Vieses AF dados")
dir()
occ <- fread(file = "Richness_XYSPP.txt", header = TRUE, fill=TRUE) %>% 
    dplyr::select(species = species,decimalLongitude = longitude,decimalLatitude = latitude) %>% 
    convert(dbl(decimalLongitude, decimalLatitude)) %>% drop_na(decimalLongitude)
save
occ %>% write.table("occ_Richness_XYSPP.txt",append = FALSE, sep = "\t", dec = ".", row.names = FALSE, col.names = TRUE,  quote = FALSE, fileEncoding = "UTF-8") #salvando
occ <- read.table("occ_Richness_XYSPP.txt", h = T, encoding = "UTF-8", sep= "\t", as.is = T)
occ <- fread("occ_Richness_XYSPP.txt", header = TRUE, fill = TRUE) %>%
    mutate(species = ifelse(species == "" | is.na(species), "non_refined", species))

occ_sf <- st_as_sf(
    occ,
    coords = c("decimalLongitude", "decimalLatitude"),
    crs = 4326,  # WGS84
    remove = FALSE  # mantém as colunas originais!
)

setwd("C:/Users/Motion/Desktop/Vieses AF dados")
dir()
shape <- st_read("shape_cortado_raster_0_25.gpkg") #abrir sua área estudo aqui (JAPA) ela tem de ser um shape quadriculado em 0.25 (qualquer duvida me pergunta)
shape <- st_transform(shape, 4326)
idx <- st_intersects(occ_sf, shape, sparse = FALSE)
occ_inside <- occ_sf[idx, ]
occ_inside_df <- as.data.frame(occ_inside)
occ <- occ_inside_df %>% as_tibble()


glimpse(occ)

setwd("C:/Users/Motion/Desktop/Vieses AF dados")
pasta <- "C:/Users/Motion/Desktop/Vieses AF dados"
dir()
dir()
arquivos_shp <- list.files(pasta, pattern = "\\.gpkg$", full.names = TRUE)
lista_shapes <- lapply(arquivos_shp, function(x) st_read(x))


plot(lista_shapes[[1]])

setwd("C:/Users/Motion/Documents/Gaz")
assentamento <- st_read("assentamento.shp")

setwd("C:/Users/Motion/Documents/Gaz")
populatedplaces <- st_read("populatedplaces.shp")#2

setwd("C:/Users/Motion/Documents/Gaz")
quilombo <- st_read("quilombo.shp")

setwd("C:/Users/Motion/Documents/Gaz")
rodestadual <- st_read("rodestadual.shp")

setwd("C:/Users/Motion/Documents/Gaz")
rodfederal <- st_read("rodfederal.shp")

setwd("C:/Users/Motion/Documents/Gaz")
terraindigenas <- st_read("terraindigenas.shp")

setwd("C:/Users/Motion/Documents/Gaz")
uc <- st_read("UC.shp")

setwd("C:/Users/Motion/Documents/Gaz")
ferrovia <- st_read("ferrovia.shp")

setwd("C:/Users/Motion/Documents/Gaz")
urbanareas <- st_read("urbanareas.shp")

data = list(ASS = assentamento, PP = populatedplaces, QU = quilombo, RES = rodestadual, REF = rodfederal, TI = terraindigenas, UC = uc, FE = ferrovia, UA = urbanareas) #colocar os outros aqui tbm, mudando as siglas pra letra maiuscula

data <- lapply(data, function(obj) {
    if (inherits(obj, "sf")) {
        return(terra::vect(obj)) 
    } else {
        return(obj)
    }
})


resolucoes <- c(0.25) #coloca 0.25
rm(resultados)
gc()
setwd("C:/Users/Motion/Desktop/Vieses AF dados")
dir()
shape <- st_read("shape_cortado_raster_0_25.gpkg")
plot(shape)
resultados <- map2(resolucoes, lista_shapes, function(res, shape) {
    print(paste("Rodando calculate_bias para resolução:", res))
    
    resultado <- calculate_bias(
        x = occ, 
        gaz = data, 
        res = res, 
        restrict_sample = shape, 
        terrestrial = TRUE
    )
    
    print(paste("Concluído para resolução:", res))
    return(resultado)
})

setwd("C:/Users/Motion/Desktop/Vieses AF dados")
getwd()

for (i in seq_along(resultados)) {
    tabela <- summary(resultados[[i]]) %>% as.data.frame()
    tabela <- cbind(factors = rownames(tabela), tabela)
    nome_arquivo <- paste0("tabela_res_", resolucoes[i], ".csv")
    write.csv(tabela, file = nome_arquivo, row.names = FALSE)
}



for (i in seq_along(resultados)) {
    
    proj <- project_bias(resultados[[i]])
    
    map_plot <- map_bias(proj, type = "log_sampling_rate")
    
    nome_mapa <- paste0("map_bias_res_", resolucoes[i], ".png")
    
    ggsave(nome_mapa, map_plot, width = 8, height = 6, dpi = 600)
    
    # rasters de intensidade de viés
    Primeiro_mais_forte  <- proj[[1]]
    Segundo_mais_forte   <- proj[[2]]
    Terceiro_mais_forte  <- proj[[3]]
    Quarto_mais_forte    <- proj[[4]]
    
    writeRaster(Primeiro_mais_forte,
                filename = paste0("raster_bias1_res_", resolucoes[i], ".tif"),
                overwrite = TRUE)
    
    writeRaster(Segundo_mais_forte,
                filename = paste0("raster_bias2_res_", resolucoes[i], ".tif"),
                overwrite = TRUE)
    
    writeRaster(Terceiro_mais_forte,
                filename = paste0("raster_bias3_res_", resolucoes[i], ".tif"),
                overwrite = TRUE)
    
    writeRaster(Quarto_mais_forte,
                filename = paste0("raster_bias4_res_", resolucoes[i], ".tif"),
                overwrite = TRUE)
    
}


plot_sampbias_part1 <- function(x, ...) {
    library(viridis)
    plo1 <- x$bias_estimate %>%
        pivot_longer(cols = contains("w_"), names_to = "bias",
                     values_to = "posterior_estimate") %>%
        mutate(bias = gsub("w_", "", .data$bias)) %>%
        mutate(bias = fct_reorder(.data$bias, .data$posterior_estimate, .fun = median, .desc = FALSE))
    p1 <- ggplot(plo1) +
        geom_boxplot(aes(x = .data$bias, y = .data$posterior_estimate, fill = .data$bias)) +
        scale_fill_viridis(discrete = TRUE) +
        xlab("Biasing factor") +
        ylab("Posterior weight") +
        coord_flip() +
        theme_bw() +
        theme(panel.grid.minor.x = element_blank(),
              panel.grid.major.y = element_blank(),
              legend.position = "none")
    print(p1)
}

plot_sampbias_part2 <- function(x, ...) {
    library(viridis)
    plo1 <- x$bias_estimate %>%
        pivot_longer(cols = contains("w_"), names_to = "bias",
                     values_to = "posterior_estimate") %>%
        mutate(bias = gsub("w_", "", .data$bias)) %>%
        mutate(bias = fct_reorder(.data$bias, .data$posterior_estimate, .fun = median, .desc = FALSE))
    plo2_w <- colMeans(x$bias_estimate)
    plo2_dist <- seq(1, 1000, length.out = 1000)
    plo2 <- data.frame(dist = plo2_dist,
                       rate = plo2_w[3] * exp(-plo2_w[4:(length(plo2_w) - 1)] * plo2_dist),
                       id = names(x$bias_estimate)[-c(1:3, ncol(x$bias_estimate))]) %>%
        mutate(id = gsub("w_", "", .data$id)) %>%
        mutate(id = factor(.data$id, levels = levels(plo1$bias)))
    plo2 <- na.omit(plo2)
    p2 <- ggplot(plo2) +
        geom_point(aes(x = .data$dist, y = .data$rate, color = .data$id)) +
        scale_color_viridis(discrete = TRUE) +
        xlab("Proximity to the bias [km]") +
        ylab("Sampling rate") +
        theme_bw() +
        theme(panel.grid.minor.x = element_blank(),
              panel.grid.major.y = element_blank(),
              legend.title = element_blank(),
              legend.position = "bottom")
    print(p2) }

for (i in seq_along(resultados)) {
    # Gerando os gráficos para cada resolução
    p1 <- plot_sampbias_part1(resultados[[i]])
    p2 <- plot_sampbias_part2(resultados[[i]])
    
    # Nome dos arquivos com a resolução incluída
    nome_p1 <- paste0("plot_part1_res_", resolucoes[i], ".png")
    nome_p2 <- paste0("plot_part2_res_", resolucoes[i], ".png")
    
    # Salvando os gráficos
    ggsave(nome_p1, p1, width = 8, height = 6, dpi = 600)
    ggsave(nome_p2, p2, width = 8, height = 6, dpi = 600)
}



#############
In case of breaking code during running PLOT2
#############

plot_sampbias_part2 <- function(x, ...) {
    library(viridis)
    plo1 <- x$bias_estimate %>%
        pivot_longer(cols = contains("w_"), names_to = "bias",
                     values_to = "posterior_estimate") %>%
        mutate(bias = gsub("w_", "", .data$bias)) %>%
        mutate(bias = fct_reorder(.data$bias, .data$posterior_estimate, .fun = median, .desc = FALSE))

    plo2_w <- colMeans(x$bias_estimate)

    plo2_dist <- seq(1, 1000, length.out = 1000)

    q <- plo2_w["q"]

    w <- plo2_w[grep("^w_", names(plo2_w))]

    ids <- gsub("w_", "", names(w))

    rates <- q * exp(-outer(plo2_dist, w))

    plo2 <- data.frame(
        dist = rep(plo2_dist, times = length(w)),
        rate = as.vector(rates),
        id = rep(ids, each = length(plo2_dist))
    )

    plo2$id <- factor(plo2$id, levels = levels(plo1$bias))

    plo2 <- na.omit(plo2)

    p2 <- ggplot(plo2) +
        geom_point(aes(x = dist, y = rate, color = id)) +
        scale_color_viridis(discrete = TRUE) +
        xlab("Proximity to the bias [km]") +
        ylab("Sampling rate") +
        theme_bw() +
        theme(panel.grid.minor.x = element_blank(),
              panel.grid.major.y = element_blank(),
              legend.title = element_blank(),
              legend.position = "bottom")

    print(p2)
}
